#ifndef CTRACKMP_H
#define CTRACKMP_H
#include "CVector3d.h"
#include <vector>
#include <string>
using namespace std;
class CBeam;
class CDtlCellDesign;
class CSpaceCharge;
class CTrackMp
{
  public:
    CTrackMp(double dt=0.1);
    ~CTrackMp();
    CTrackMp(double,double dt=0.1);
    void           DtlCell(double,double,double,double,double,double,double,double,double,double,double,int);
    void           Drift(double);
    void           DTL_Fourier(string,double);
    void           DTL_Fourier_2(string,double,int);
		double         Drift_Matrix(double);
		void           Quad_Matrix(double dQL,double dQS);
    void           Quad(double,double);
    void           ThinGap(double,double,double,double,double);
    void           OutCordFile(int);
    void           SetFreq(double);
    void           SetDt(double);        //输入单位是ns   //这个函数存在没意义？？？beam里面有一个就可以了
    void           GetParticleNum();
    void           SetSfline(vector<double> *);
    void           SetBetag(double);
    void           SetTTF(double,double,double,double);
    void           SetCellDesign(CDtlCellDesign* _cDtlCell) { cDtlCell=_cDtlCell; };
    void           SetPhaseMove(double _dPhi);
    void           OutBeamRmsSize();
    void           OutTwiss();
    void           getBeam();
    void           SetNSC(int _NSC=1)  {NSC=_NSC;};//表示每几步进行一次SC的计算
    void           Read_FieldMap_EMCyl(string); 
    void           QuadForFieldMap(string);//临时用来存储quad信息 
    void           DtlMap(string);
    void           MapForTW(string);
    void           SParticleTrack();
    void           SearchInitPhase();
    void           DTL_FieldMap(string,string,double);
    void           SetBoreRadius(double _dBoreRadius=0);
    int            CheckBeamSize();
		double         FringeQuad(double,double,double);
    void           gsl_test();  //temperary
		vector<double> getMatrixInput_Twiss();
		vector<double> getMatrixInput_Beam();
		void           testParticle();
//		void           calSig0();

  private:
    double         XI0(double);
    double         XI1(double);
    void           Interpret(double,double,double,CVector3d *,CVector3d *);
    double         InterpretQuad(double);
    void           GapPhase(double,double);//find the Tank input phase according to the synchrotron phase in the first cell.
    // Need to check whether the cell is the first one of the whole accelerator 
  public:
    CBeam         *beam;
    CSpaceCharge  *SC;
    int           iParticleNumber;
    int           nCell;
    double        dT,dTp;
    double        drMap,dzMap;      //step size
    double        drMin,drMax,dzMin,dzMax; //sacle
    int           NrMap,NzMap;             //step numbers
    double        dBoreRadius;
    double        dt;
		double        dZedge;
	//	double        dSig0r,dSig0z; //phase advance
  private:
    int           NSC;
    CDtlCellDesign *cDtlCell;
    double        dCellLength,dLeft,dRight,dGap,dLQuad1,dLQuad2,dSQuad1,dSQuad2;
    double        dE0,dPhiS,dTTF,dTTFp,dSTF,dSTFp;  //transit time factor T and S
    double        dAperture;
    double        dBeta,dBetac,dGammac,dEnergyc,dPhic;
    double        dEnergyS,dEnergyi,dEnergyES,dGammai,dGammaES,dBetai,dBetaES,dBetag,dBetaS,dGammaS;
    double        dWaveLength,dFreq;
    double        dX,dY,dZ;
    int           iFlagThinGap;                
    double        dPhaseChange;
    vector<vector<double> > vvdBeamRms;
    vector<vector<double> > vvdTwissRms;
    vector<vector<vector<double> > > vvvdFMap;
    vector<vector<double>  > vvdQuadMap;
    vector<vector<CVector3d> > vvdEMap;
    vector<vector<CVector3d> > vvdBMap;
    int          iOld;
};
#endif
